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Abstract 

This work presents the application of density-based topology optimisation 
to the design of three-dimensional heat sinks cooled by natural convection. 
The governing equations are the steady-state incompressible Navier-Stokes 
equations coupled to the thermal convection-diffusion equation through the 
Bousinessq approximation. The fully coupled non-linear multiphysics system 
is solved using stabilised trilinear equal-order finite elements in a parallel 
framework allowing for the optimisation of large scale problems with order of 
40-330 million state degrees of freedom. The flow is assumed to be laminar 
and several optimised designs are presented for Grashof numbers between 
10^ and 10®. Interestingly, it is observed that the number of branches in 
the optimised design increases with increasing Grashof numbers, which is 
opposite to two-dimensional optimised designs. 

Keywords: topology optimisation, heat sink design, natural convection, 
large scale, multiphysics optimisation 


1. Introduction 

Natural convection is the phenomena where density-gradients due to tem¬ 
perature differences cause a fluid to move. Natural convection is therefore 
a natural way to passively cool a hot object, such as electronic components, 
light-emitting diode lamps or materials in food processing. 
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Topology optimisation is a material distribution method [l| used to op¬ 
timise the layout of a structure in order to minimise a given performance 
measure subject to design constraints and a physical model. In order to take 
convective heat transfer, to an ambient fluid, into account in the design pro¬ 
cess of density-based methods, a common extension is to introduce some form 
of interpolation of the convection boundaries, see e.g. @,0,1^. More recently, 
Dede et al. used these simplihed models to design and manufacture heat 
sinks subject to jet impingement cooling. However, it is hard to justify the 
application of a predetermined and constant convection coefficient, because 
topology optimisation often leads to unanticipated designs, closed cavities 
and locally varying velocities. During the optimisation process, the design 
also changes significantly and the interaction with the ambient fluid changes. 
Therefore, to ensure physically correct capturing of the aspects of convective 
heat transfer, the full conjugate heat transfer problem must be solved. 

Topology optimisation for fluid systems began with the treatment of 
Stokes flow in the seminal article by Borrvall and Petersson and has 
since been applied to Navier-Stokes j^, as well as scalar transport prob¬ 
lems amongst others. The authors have previously presented a density- 
based topology optimisation approach for two-dimensional natural convec¬ 
tion problems j^. Recently, Coffin and Maute presented a level-set method 
for steady-state and transient natural convection problems using X-FEM 
Tol l. Interested readers are referred to for further references on topology 


optimisation in fluid dynamics and heat transfer. 

Throughout this article, the flows are assumed to be steady and laminar. 
The fluid is assumed to be incompressible, but buoyancy effects are taken 
into account through the Boussinesq approximation, which introduces vari¬ 
ations in the fluid density due to temperature gradients. The inclusion of a 
Brinkman friction term facilitates the topology optimisation of the fluid flow. 

The scope of this article is primarily to present large scale three-dimensional 
results using the formulation presented in j^, as well as the computational 
issues arising from solving the non-linear equation system and the subsequent 
linear systems. Thus, only a brief overview of the underlying finite element 
and topology optimisation formulations are given and the reader is referred 
to for further information. 

In recent years an increasing body of work has been published on efficient 
large scale topology OTtimisation. These works cover the use of high-level 


scripting languages [ll|, ll2j, multiscale/-resolution approaches [13|, [l^ and 


parallel programming using the message parsing interface and C/Fortran 
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lil, 0, H O . To facilitate the solution to truly large scale conjugate heat 
transfer problems, the implementation in this paper is done using PETSc 


18 


19| and the framework for topology optimisation presented in 

The layout of the paper is as follows: Section [2] presents the governing 
equations; Section |3] presents the topology optimisation problem; Section 
m briefly discusses the hnite element formulation; Section E] discusses the 
numerical implementation details; Section E] presents scalability results for 
the parallel framework and optimised designs for a test problem; Section [ 7 ] 
hnishes with a discussion and conclusion. 


2. Governing equations 

The dimensionless form of the governing equations have been derived 
based on the Navier-Stokes and convection-diffusion equations under the as¬ 
sumption of constant fluid properties, incompressible, steady flow and ne¬ 
glecting viscous dissipation. Furthermore, the Boussinesq approximation has 
been introduced to take density-variations due to temperature-differences 
into account. A domain is decomposed into two subdomains, hi = hlj U 
where Qf is the fluid domain and is the solid domain. In order to facilitate 
the topology optimisation of conjugate natural convective heat transfer be¬ 
tween a solid and a surrounding fluid, the equations are posed in the unihed 
domain, hi, and the subdomain behaviour is achieved through the control of 
coefficients. The following dimensionless composite equations are the result. 
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where Ui is the velocity field, p is the pressure field, T is the temperature field, 
Xi denotes the spatial coordinates, ef is the unit vector in the gravitational 
direction, a{x.) is the spatially-varying effective impermeability, A'(x) is the 
spatially-varying effective thermal conductivity, s(x) is the spatially-varying 
volumetric heat source term, Pr is the Prandtl number, and Gr is the Grashof 
number. 
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The effective thermal conductivity, K{x.), is defined as: 


K{^) 


j 1 if X G 

li ifxGf^. 


( 4 ) 


where Ck = ^ is the ratio between the fluid thermal conductivity, kf, and the 
solid thermal conductivity, kg. Theoretically, the effective impermeability, 
q;(x), is defined as: 


j 0 if X G ff/ 

1 CX3 if X G ffs 


( 5 ) 


in order to ensure zero velocities inside the solid domain. However, numeri¬ 
cally this requirement must be relaxed as will be described in section [Hi The 
volumetric heat source term is defined as being active within a predefined 
subdomain of the solid domain, uj C 


s(x) 


f 0 if X ^ a; 
[so if X G a; 


( 6 ) 


where sq = dimensionless volumetric heat generation, q is the 

dimensional volumetric heat generation, AT is the reference temperature 
difference and L is the reference length scale. 

The Prandtl number is defined as: 

Pr = ^ (7) 


where v is the kinematic viscosity, or momentum diffusivity, and T is the 
thermal diffusivity. It thus describes the relative spreading of viscous and 
thermal effects. The Grashof number is defined as: 


Gr 


g(5 AT 


( 8 ) 


where g is the acceleration due to gravity and (3 is the volumetric coefficient of 
thermal expansion. It describes the ratio between the buoyancy and viscous 
forces in the fluid. The Grashof number is therefore used to describe to what 
extent the flow is dominated by natural convection or diffusion. For low 
Gr the flow is dominated by viscous diffusion and for high Gr the flow is 
dominated by natural convection. The problems in this article are assumed 
to have large enough buoyancy present to exhibit natural convective effects, 
but small enough Gr numbers to exhibit laminar fluid motion. 
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3. Optimisation formulation 


3.1. Interpolation functions 

In order to perform topology optimisation, a continuons design field, 7 (x), 
varying between 0 and 1 is introdnced. Pure fluid is represented by 7 (x) = 1 
and solid by 7 (x) = 0. For intermediate values between 0 and 1, the effective 
conductivity is interpolated as follows: 


K{3) 


7(C'fc(l + g/) - 1) + 1 

<^^(1 + g/7) 


( 9 ) 


and likewise the effective impermeability is interpolated using: 

a (7) = a /, ^ (10) 

1 + g«7 

The interpolation functions ensure that the end points dehned in (jl]) and (jSj), 
respectively, are satished. The effective impermeability is bounded to a in the 
solid regions and this upper bound should be chosen large enough to provide 
vanishing velocities, but small enough to ensure numerical stability. The 
convexity factors, qf and Qa, are used to control the material properties for 
intermediate design values in order to promote well-dehned designs without 
intermediate design held values. 


3.2. Optimisation problem 

The topology optimisation problem is dehned as: 


minimise: / ( 7 , T) = 
-yev 


sU)TdV 


subject to: g{')) = 1 — j dV <Vf dV 

Jfld 

'^(7,u,p,T) = 0 
0 < 7 (x) <1 Vx e f2d 


( 11 ) 


where 7 is the design variable held, V is the design space, / is the ther¬ 
mal compliance objective functional, g is the volume constraint functional, 
^( 7 , u,p, T) is the residual arising from the stabilised weak form of the gov¬ 
erning equations, and ^ H is the design domain. 

The design held is regularised using a PDE-based (partial diherential 


equation) density hlter [20 


the method of moving asymptotes (MMA) [2l| 


18l | and the optimisation problem is solved using 

13. 
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3.3. Continuation scheme 


A continuation scheme is performed on various parameters in order to 
stabilise the optimisation process and to improve the optimisation results. 
It is the experience of the authors that the provided continuation scheme 
yields better results than starting with the end values, although this cannot 


generally be proven [22|, [23 


The chosen continuation strategy consists of hve steps: 


Qf G {0.881, 8.81, 88.1, 881, 881} (12a) 

ga e{8,8,8,98,998} (12b) 

a e{l0^10^10^10^10^} (12c) 

The sequence is chosen in order to alleviate premature convergence to poor 
local minimum. The value of g/ is slowly increased to penalise intermedi¬ 
ate design held values with respect to conductivity. The maximum effective 
permeability, a, is set relatively low during the hrst three steps, as this en¬ 
sures better scaled sensitivities and more stable behaviour. Over the last two 
steps, a is increased by two orders of magnitude in order to further decrease 
the velocity magnitudes in the solid regions. The particular values of g/ are 
chosen by empirical inspection such as to ensure the approximate collocation 
of the huid and thermal boundaries. A more comprehensive investigation 
into the modelling accuracy of the boundary layers is outside the scope of 
this paper and will be investigated in future work. 

It is important to note that the optimisation problem is by no means 
convex and any optimised design will at best be a local minimum. The ob¬ 
tained design will always depend on the initial design as well as the continu¬ 
ation strategy. However, in the authors experience, the chosen continuation 
strategy gives a good balance between convergence speed, hnal design per¬ 
formance and physicality of the modelling. The effects of the steps of the 
current continuation strategy on the design distribution will be discussed in 
section 16.41 


4. Finite element formulation 

The governing equations are discretised using stabilised trilinear hexahe- 
dral hnite elements. The design held is discretised using elementwise con¬ 
stant variables, in turn rendering the effective thermal conductivity and the 
effective impermeability to be elementwise constant. The monolithic hnite 
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element discretisation of the problem ensures continuity of the temperature 
held, as well as the huxes across huid-solid interfaces. The particularities of 
the implemented hnite element formulation are as detailed in j^. However, 
simpler stabilisation parameters have been used in order to allow for consis¬ 
tent sensitivities to be employed, see Appendix B The Jacobian matrix is 
now fully consistent in that variations of the stabilisation parameters with 
respec_^ to design and state helds are included, in contrast to the original 
work 


5. Numerical implementation 


The discretised FEM equation is implemented in PETSc [19|] based on the 
topology optimisation framework presented in jl^. The PETSc framework 
is used due to its parallel scalability, the availability of both linear and non¬ 
linear solvers, preconditioners and structured mesh handling possibilities. All 
components described in the following are readily available within the PETSc 
framework. 


5.1. Solving the non-linear system 

The non-linear system of equations is solved using a damped Newton 
method. The damping coefficient is updated using a polynomial L^-norm £t, 
where the coefficient is chosen as the minimiser of the polynomial £t. The 
polynomial £t is built using the L^-norm of the residual vector at the current 
solution point, at 50% of the Newton step and at 100% of the Newton step. 
This residual-based update scheme combined with a good initial vector (the 
solution from the previous design iteration) has been observed to be very 
robust throughout the optimisation process for the moderately non-linear 
problems treated. To further increase the robustness of the non-linear solver, 
if the Newton solver fails from the supplied initial vector, a ramping scheme 
for the heat generation magnitude is applied in order to recover. Throughout, 
the convergence criteria for the Newton solver is set to a reduction in the L^- 
norm of the residual of 10“^ relative to the initial residual. 


5.2. Solving the linear systems 

Due to the large scale and three-dimensional nature of the treated prob¬ 
lems, the arising linearised systems of equations are by far the most time 
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consuming part of the Newton scheme. Therefore, to make large scale prob¬ 
lems tractable the (unsymmetric) linear systems are solved using a fully 
parallelised iterative Krylov subspace solver. 

Constructing an iterative solver that is both independent of problem set¬ 
tings and possesses both parallel and numerical scalability, is intricate and 
beyond the scope of this work, where focus is on the optimisation. However, 
in order to facilitate the solution of large scale optimisation problems, an 
efficient solver is required. To this end the authors use the readily available 
core components in PETSc to construct a solver with focus on simplicity, as 
well as reduced wall clock time. This is quite easy to obtain for solvers and 
preconditioners that rely heavily on matrix-vector multiplications. To this 
purpose the flexible GMRES (F-GMRES) method is used as the linear solver 
combined with a geometry-based Galerkin-projection multigrid (GMG) pre¬ 
conditioner. Such a solver depends highly on the quality of the smoother 
to guarantee fast convergence [^. The authors have found that a simple 
Jacobi-preconditioned GMRES provides a reasonable choice of smoothei0. 
The convergence criterion for the Krylov solver is set to 10“® relative to the 
initial residual. 

6. Results 

6.1. Problem setup 

The considered optimisation problem is an (academic) example of a heat 
sink in a closed cubic cavity. Figure [1] shows the problem setup with di¬ 
mensions and boundary conditions. The heat source (black in hgure), exem¬ 
plifying an electronics chip, is placed in the mid-bottom of the cavity and 
modelled using a small block of solid material with volumetric heat genera¬ 
tion, So = lO'^. The design domain (gray in hgure) is placed on top of the 
heat source in order to allow the cooling huid to pass underneath it, as well 
as to allow room for wires etc. The vertical and top outer walls of the cavity 
are assumed to be kept at a constant cold temperature, T = 0, while the 
bottom is insulated. The height of the cavity has been used as the reference 
length scale, L. Thus, the cavity dimensions are 1x1x1, the design do¬ 
main dimensions are 0.75 x 0.75 x 0.75 and the heat source dimensions are 


^Due to the choice of a Krylov smoother, the multigrid preconditioner will vary slightly 
with input and thus require a flexible outer Krylov method. 



U, = OAT = 0 





Ui =0 AT = 0 


n 


Wi = 0 A —K^Ui = 0 


(a) Dimensions (b) Boundary condi¬ 
tions 


Figure 1: Illustration of the problem setup for the heat sink in a close cubic 
cavity. The heat source is black and the design domain is gray. 


0.1 X 0.05 X 0.1. A discussion on the dehnition of the temperature difference 
in the Grashof number can be found in Appendix A Initial investigations 
showed that due to the symmetry of the domain and boundary conditions, 
the design and state solutions remained quarter symmetric throughout the 
optimisation. Thus, the computational domain is limited to a quarter of the 
domain with symmetry boundary conditions. The volume fraction is set to 
5%, i.e. Vf = 0.05, for all examples. 


6.2. Parallel performance 

To demonstrate the parallel performance of the state solver, the model 
optimisation problem described in section [6] is solved on a hxed mesh for dif¬ 
ferent numbers of processes. All computations in this paper were performed 
on a cluster, where each node is equipped with two Intel Xeon e5-2680v2 10- 
core 2.8GHz processors. The results shown in Tables [U and [2] are averaged 
over 250 design cycles and show the performance for Gr = 10^ and Gr = 10®, 
respectively. The data shows that the proposed solver scales almost linearly 
in terms of speed up, and more importantly that the performance is only 
slightly affected by the Grashof number. 

In order to quantify the degree of numerical scalability, a second study is 
performed in which the mesh resolution is varied. The study is conducted for 
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Processes 

time [s 

scaling 

160 

53.2 

1.00 

320 

28.9 

0.54 

640 

14.1 

0.26 


Table 1: Average time taken per state solve over 250 design iterations for 
Gr = 10^ at a mesh resolution of 80 x 160 x 80. 


Processes 

time [s 

scaling 

160 

62.6 

1.00 

320 

31.9 

0.51 

640 

16.5 

0.26 


Table 2: Average time taken per state solve over 250 design iterations for 
Gr = 10® at a mesh resolution of 80 x 160 x 80. 


Mesh size 

Gr = 10® 

Gr = 10® 

80 X 160 X 80 

7.5 

5.6 

160 X 320 X 160 

10.1 

7.7 

320 X 640 X 320 

18.4 

15.6 


Table 3: Average iterations for the linear solver per state solve over entire 
design process for Gr = 10® and Gr = 10® at varying mesh resolutions. 
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both low and high Grashof numbers and the average F-GMRES iterations 
are collected in Table |31 The total number of design iterations averaged 
over was 250, 500 and 1000, respectively, for the three mesh resolutions. 
The data clearly shows that the computational complexity increases with 
problem size, and thus that the solver is not numerically scalable. However, 
since the growth in numerical effort, i.e. the number of F-GMRES iterations, 
is moderate we conclude that the proposed solver is indeed applicable for 
solving large scale natural convection problems. 


6.3. Varying Grashof number 

The problem is investigated for varying Gr under constant volumetric 
heat generation, sq = 10^) Prandtl number, Pr = 1, and thermal conduc¬ 
tivity ratio, Ck = 10“^. The purpose of the present study is not to provide 
a detailed physical example. It is rather to provide phenomenological in¬ 
sight into the effect of changing the governing parameter of the fluid-thermal 
coupling, namely the Grashof number. However, physical interpretations of 
tuning the Gr-number, while keeping the dimensionless volumetric heat gen¬ 
eration and the Pr-number constant, could be e.g. equivalent to tuning the 
gravitational strength (going from microgravity towards full gravity) or the 
dimensional volumetric heat generation. It is important to note that when 
interpreting the results, the dimensional temperature scale will differ for the 
two interpretations. While tuning the gravitational strength, the tempera¬ 
ture scale remains the same; by varying the magnitude of the volumetric heat 
generation, the temperature scale varies accordinglj{§. 

The computational mesh is 160 x 320 x 160 elements yielding a total 
of 8,192, 000 elements and 41, 603, 205 degrees of freedom for the quarter 
domain. The design domain consists of 3,456,000 elements and the hlter 
radius is set to 2.5 times the element size. 

Figure |2] shows the optimised designs for varying Gr-number with super¬ 
imposed slices of the corresponding temperature helds. Due to the use of 
density hltering, the interface between solid and fluid regions for the opti¬ 
mized designs are not exactly described but consists of a thin layer of inter¬ 
mediate design held values. The optimised design distributions are shown 


^The dimensionless volumetric heat generation is given by sq = Requiring that 

So remains constant, means that the dimensional volumetric heat generation, q, and the 
reference temperature difference, AT, must vary in unison, i.e. = const.. An increase 
in the Gr-number is thus achieved through the increase of AT. 
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(a) Gr = 103 (b) Gr = 10^ 


Temperature 

.860435 

; 0.8 


0.75 

0.5 


jO.25 

O' 



0 


,V! 


(c) Gr = 10^ 


(d) Gr = 10® 


Figure 2: Optimised designs for varying Gr-number at a mesh resolution of 
160 X 320 X 160. 



Temperature 

.1638 
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Gr 

Primary 

Secondary 

Surface area 

10 ® 

12 

0 

0.887 

10 ^ 

8 

16 

0.853 

10 ® 

8 

28 

0.834 

10 ® 

8 

48 

0.846 


Table 4: The number of primary and secondary branches, as well as the 
surface area for the optimised designs of figure [ 2 l 


thresholded at 7 = 0.05, which is the approximate location of the computa¬ 
tional interface. The general characteristics of all the designs are similar, i.e. 
all designs are “thermal trees” with conductive branches moving the heat 
away from the heat source. However, it can clearly be seen that the design 
changes considerably with increasing convection-dominance (increasing Gr). 
For increasing Gr-number the conductive branches contract, resulting in a 
smaller spatial extent of the overall heat sink. This intuitively makes sense as 
the problem goes from one of conduct ion/diffusion at Gr = 10^ to convection 
at Gr = 10®. When diffusion dominates, the goal for the branches essentially 
becomes to conduct the heat directly towards the cold walls. As convection 
begins to matter, the fluid movement aids the transfer of heat away from the 
heat sink and the branches do not need to be as long. Instead, the design 
forms higher vertical interfaces in order to increase surface area perpendicu¬ 
lar to the flow direction and thus better transfer of heat to the fluid. At the 
same time the complexity of the designs increases as the importance of con¬ 
vection increases. This can be seen by studying the number of primary and 
secondary branches. Primary branches are dehned as those connected to the 
heat source directly and secondary branches are those connected to primary 
branches. Table 0] shows the number of primary and secondary branches of 
the optimised designs shown in figure 01 The number of primary branches is 
largest for the diffusion-dominated case, but more or less constant thereafter. 
However, it can be seen that the number of secondary branches significantly 
increases as the Gr-number is increased. Table 0] also shows that the total 
surface areals is decreasing for the three lower Gr-numbers and then increases 
slightly. 


^The surface area is computed using an isosurface at the selected threshold design field 
value. 
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Figure 3: Temperature distribution in optimised designs for Gr = 10^ and 
Gr = 10® at a mesh resolution of 160 x 320 x 160 - view from below. 


It is interesting to note, that the trend of increasing complexity with 
increasing Gr-number is the reverse of what was observed for two-dimensional 
problems . There, the complexity of the design decreased as the Gr-number 
increased. This difference is likely due to the fact that going into three- 
dimensions allows the fluid to move around and through the design making 
it more a question of “topology”, in contrast to in two-dimensions where it is a 
question of surface shape. Physically, in two-dimensions additional branches 
disturb the flow and thus the heat transfer; in three-dimensions, vertical 
branches improve the heat transfer through an increased vertical surface 
area. Figure [3] shows the optimised designs for Gr = 10^ and Gr = 10® 
from below. The radial extent of the designs are emphasised from these 
views. It is also seen that the branches for Gr = 10® are positioned to ensure 
that the structure is open from below, with the branches forming vertical 
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Optimisation Gr 


Analysis Gr 

10® 

10^ 

10® 

10® 

10® 

2.064962 

2.066856 

2.241125 

2.362470 

10^ 

1.933460 

1.880946 

1.994538 

2.111726 

10® 

1.488278 

1.450276 

1.404511 

1.439089 

10® 

1.134963 

1.121979 

1.061452 

1.025340 


Table 5: Crosscheck objective function values for the designs shown in figure 

H 


Gr Time Non-linear: Linear: 

avg. (max) contin - avg. (max) 


10® 

9:56 

1.9 (2) 

7.6,8.4,7.8,8.1,18.4 - 

10.1 (25) 

10^ 

10:25 

2.0 (3) 

8.3,8.6,8.3,8.6,22.7- 

11.3 (29) 

10® 

10:28 

2.1 (10) 

8.4,8.6,8.7,8.2,15.7 

- 9.9 (34) 

10® 

10:35 

2.1 (7) 

7.3,7.4,7.5,8.0,8.4 - 

7.7 (14) 


Table 6: Computational time, average non-linear iterations and linear itera¬ 
tions for the optimised designs of figure [2l 


walls as discussed above. 

Table [5] shows a crosscheck of the objective functions for the optimised 
designs. It can be seen that all designs optimised for certain flow conditions 
outperform the other designs for the specified Gr-number. 

The optimisations was run for 500 design iterations for each optimised 
design. Table |6] shows the computational time, average non-linear iterations 
and linear iterations for the optimised designs of figure |2] using 1280 cores. 
It can be seen that the computational time only weakly increases as the Gr- 
number is increased and remains close to 1.2 minutes per design iteration 
on average. It is interesting to note, that Gr = 10® seems easier to solve 
than the others as it exhibits lower average number of linear iterations than 
all others, as well as a lower maximum number of non-linear iterations than 
Gr = 10®. Furthermore, it can be seen that due to the Newton method 
starting from a good initial vector, only 2 non-linear iterations are needed 
for most of the design iterations independent of Gr-number. 
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6.4- High resolution design 

The problem is now investigated with a computational mesh of 320 x 
640 X 320 elements yielding a total of 65, 536, 000 elements and 330, 246,405 
degrees of freedom for the quarter domain. The design domain consists of 
27, 648, 000 elements and the filter radius is set to 2.5 times the element size, 
i.e. in absolute measures half the size of before. The optimisation is run for 
1000 design iterations and the computational time was 107 and 108 hours, 
respectively, for Gr = 10^ and Gr = 10®, using 2560 cores. This yields an 
average of 6.4 and 6.5 minutes per design iteration, respectively. 

Figure H] shows the optimised design for Gr = 10® with the hne mesh 
resolution and small length scale at various steps of the continuation strat¬ 
egy. The complexity of the design can be seen to be signihcantly higher than 
for the design with a larger length scale, figure |2dl The subfigures are the 
hnal iterations of the 1st, 2nd, 3rd and 5th (final) continuation steps. It 
can be seen that the complexity of the design decreases during the optimi¬ 
sation process once the overall topology has been settled (iteration 400 and 
onwards). This is due to the harder and harder penalisation of intermediate 
design held values, with respect to conductivity, which are present at the 
interface between solid and huid. Therefore, smaller features are progres¬ 
sively removed as the surface area is more heavily penalised. The reason 
for going to such high penalisation of the conductivity is to ensure the ap¬ 
proximate collocation of the huid and thermal boundaries. However, if one 
starts directly with these physically-relevant parameters, particularly poor 
local minima have been observed. 

Figure |5] shows the hnal optimised designs for Gr = 10^ and Gr = 10® 
with the hne mesh resolution and small length scale. The complexities of 
both designs can be seen to be signihcantly higher than the previous, hgures 
Ealand IMl due to the smaller length scale. The increase in complexity with 
increasing Grashof number is not as apparent at this resolution and length 
scale, however, it is still argued that the number of primary and secondary 
members follows the trend from the previous study, section 16.31 Also, the 
shift from long conducting branches to shorter members, with signihcant 
vertical surface area, is still obvious. 

7. Discussion and conclnsion 

In this paper we have applied topology optimisation to the design of 
three-dimensional heat sinks using a fully coupled non-linear thermohuidic 
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(c) Iteration 600 (d) Iteration 1000 


Figure 4: Optimised designs for Gr = 10® at a mesh resolution of 320 x 
640 X 320. Please note that the freely hanging material in (a) is due to only 
elements below the threshold, 7 = 0.05, are shown - the design is connected 
by intermediate design held values. 
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(b) Gr = 10® 


(a) Gr = 10^ 


Figure 5: Optimised designs for Gr = 10^ and Gr = 10® at a mesh resolution 
of 320 X 640 X 320. The outline of the outer walls are shown in black. 


model. In contrast to previous works, that considered simplihed convection 
models, the presented methodology is able to recover interesting physical 
effects and insights, and avoids problems with the formation of non-physical 
internal cavities, length-scale effects and artihcial convection assumptions. 
The implementation of the code in a PETSc framework suitable for large 
scale parallel computations allows for running examples with more than 300 
million degrees of freedom and almost 30 million design variables on regular 
grids. 

The example considered in the paper was primarily of academic nature. 
Nevertheless, some interesting insight is obtained, showing that optimised 
structures go from exhibiting simple branches that conduct heat towards the 
cold outer boundaries for diffusion-dominated problems, towards complex 
and compact multi-branched structures that maximize the convection heat 
transfer for higher Grashof numbers. 

Current and future work includes applications to real life problems (see 
preliminary work in 2^), irregular meshes, multiple orientations, as well 
as the extension to transient problems and detailed investigations into the 
modelling accuracy of the boundary layer. 


18 










Computational Gr Mesh Figure Actual Gr 


10^ 

160 

X 

320 

X 

160 

[2al 

1.69 

X 

10^ 

10^ 

160 

X 

320 

X 

160 

[2b] 

1.55 

X 

10^ 

10^ 

160 

X 

320 

X 

160 

M 

1.16 

X 

10® 

10® 

160 

X 

320 

X 

160 

M 

8.60 

X 

10® 

10^ 

320 

X 

640 

X 

320 

Ibal 

1.42 

X 

10^ 

10® 

320 

X 

640 

X 

320 

[5b] 

7.16 

X 

10® 


Table A.7: Computational and actual Gr-numbers for the designs presented 
in this paper. 
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Appendix A. Computational versus actual Grashof number 

The Grashof numbers used above are all based on an a priori dehned ref¬ 
erence temperature difference. However, the actual temperature difference 
observed between the heat source and the walls are not known a priori due 
to the fact that the problem only has a single known temperature (Dirichlet 
boundary condition) and a volumetric heat source. This is why the maxi¬ 
mum temperature observed for the designs, see hgure [2], is not equal to 1. 
Therefore, after the optimisation, one can dehne an a posteriori Grashof 
number based on the actual computed temperature difference. The a priori 
version can be termed the computational Grashof number, the one that goes 
into the dimensionless governing equations; and the a posteriori version can 
be called the actual, or optimised, Grashof number for the given optimised 
design. The actual Grashof number can be useful for experimental studies, as 
well as for future comparisons. The actual Grashof numbers for the designs 
shown in this paper are listed in Table I A. 71 The values for the hne meshes 
are lower due to the lower thermal compliance, equivalent to the temperature 
of the heat source, allowed by the smaller design length scale. 
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Appendix B. Stabilisation parameters 


The stabilised weak form of the equations are given in [9| . The equations 
are stabilised using the pressure-stabilising Petrov-Galerkin (PSPG) j^, 27 
and streamline-upwind Petrov-Galerkin (SUPG) methods [28|, for more in¬ 
formation please see j^. 

The stabilisation parameters for SUPG and PSPG are dehned as one and 
the same using the following approximate min-function: 


TSU = Tps = T 




1 

2 


The limit factors are given by: 


U 


7-2 


TS 


4he 


Ug 2 

APr 

1 

Og 


(B.l) 


(B.2a) 

(B.2b) 

(B.2c) 


where he is a characteristic element size (for cubes the element edge length) 
and Ug is the element vector of velocity degrees of freedom. The first limit 
factor, Ti, has been simplified based on evaluation at element centroids un¬ 
der the assumption of a single Gauss-point, yielding a constant stabilisation 
factor within each element. 

In order to define a consistent Jacobian matrix, and thus a consistent 
adjoint problem, the derivatives of the stabilisation factors are needed with 
respect to the velocity held. This can be found to be: 


dr 

dvie 



(B.3) 


Furthermore, to dehne consistent design sensitivities, the derivatives of 
the stabilisation factors with respect to the design held is needed. This can 
be found to be: 


dr T 

dje n 




(B.4) 
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Including the derivatives in the dehnition of consistent adjoint sensitivities 
has been observed to provide a one to two order of magnitude improvement 
in accuracy of sensitivities with respect to a hnite difference approximation. 
Signihcant differences have not been observed in optimisation behaviour or 
in hnal designs. However, this cannot be guaranteed in general and it is 
therefore best to ensure consistency, as also highlighted by the similar issue 
of frozen turbulence 
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Similar dehnitions and derivations are carried out for the thermal SUPG 
stabilisation. 
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